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We have studied the dynamics of an equal-mass magnetized neutron-star binary within a resistive magnetohy¬ 
drodynamic (RMHD) approach in which the highly conducting stellar interior is matched to an electrovacuum 
exterior. Because our analysis is aimed at assessing the modifications introduced by resistive effects on the 
dynamics of the binary after the merger and through to collapse, we have carried out a close comparison with 
an equivalent simulation performed within the traditional ideal magnetohydrodynamic approximation. We have 
found that there are many similarities between the two evolutions but also one important difference: the survival 
time of the hypermassive neutron star increases in a RMHD simulation. This difference is due to a less efficient 
magnetic-braking mechanism in the resistive regime, in which matter can move across magnetic-field lines, 
thus reducing the outward transport of angular momentum. Both the RMHD and the ideal magnetohydrody¬ 
namic simulations carried here have been performed at higher resolutions and with a different grid structure than 
those in previous work of ours [L. Rezzolla, B. Giacomazzo, L. Baiotti, J. Granot, C. Kouveliotou, and M. A. 

Aloy, Astrophys. J. Letters 732, L6 (2011)], but confirm the formation of a low-density funnel with an ordered 
magnetic field produced by the black hole-torus system. In both regimes the magnetic field is predominantly 
toroidal in the highly conducting torus and predominantly poloidal in the nearly evacuated funnel. Reconnection 
processes or neutrino annihilation occurring in the funnel, none of which we model, could potentially increase 
the internal energy in the funnel and launch a relativistic outflow, which, however, is not produced in these 
simulations. 

PACS numbers: 04.25.Dm, 04.40.Dg, 95.30.Lz, 97.60.Jd 

I. INTRODUCTION 

With the rapid progress made in upgrading and testing a se¬ 
ries of advanced interferometric gravitational-wave detectors 
such as LIGO [ ], Virgo [2], and KAGRA [3], there are now 
great expectations that in the next five years we will witness 
the first direct detection of gravitational waves. Prime sources 
for such a detection are binary systems of compact objects, 
namely, binary systems comprising either two black holes, a 
black hole and a neutron star, or two neutron stars. The lat¬ 
ter configuration, in particular, is potentially a very interest¬ 
ing one, as it will represent the most common source, with a 
realistic expected detection rate of ^ 40yr“^ [4]. A detec¬ 
tion of gravitational waves from binary neutron stars would 
yield a wealth of information about the chirp mass, the ori¬ 
entation, and the localization of the binary but also possibly 
the mass, spin and radius of the individual stars [5, 6]. In 
turn, this information could set constraints on the equation of 
state (EOS) of the matter in their interior. Indeed, a number 
of recent investigations have revealed that it is possible to set 
serious constraints on the properties of the neutron-star struc¬ 
ture and EOS, either when using the inspiral signal only [7, 8], 
or when exploiting the rich spectral features of the postmerger 
signal [9-13]. 

At the same time, the merger of a binary system containing 
at least one neutron star represents arguably the most attrac¬ 
tive scenario to explain the complex phenomenology associ¬ 
ated with short gamma-ray bursts (SGRBs), although many 
alternatives exist [see Ref. [14] for a recent review]. While 
such a scenario was suggested already 30 years ago [15, 16], 
numerical simulations (see, e.g.. Refs. [17-21]) and new ob¬ 


servations [14, 22] have put this scenario on firmer grounds. 
In particular, what these simulations have shown is that the 
merger of a binary system of neutron stars inevitably leads 
to the formation of a metastable object, which we dub the 
binary-merger product or BMP. Depending on the total mass 
and mass ratio of the binary, the BMP can either be a supra- 
massive neutron star (SMNS), that is, a star with mass above 
the maximum mass for nonrotating configurations but 

below the maximum mass for uniformly rotating configura- 
tions Mmax, with M^ax ^ (1-15 - 1.20) [23]; a hy- 

permassive neutron star (HMNS), that is, a star more massive 
than a SMNS; or a black hole^. 

The general-relativistic hydrodynamical modelling of bi¬ 
nary neutron stars has seen very considerable progress over 
the last decade (see, e.g.. Refs. [17-19]), and it has now 
reached a rather mature state. In fact, it is presently possible 
to calculate inspiral waveforms having a phase accuracy com¬ 
parable to that of binary black-hole simulations thanks to the 
use of high-order methods with high-order convergence rates 
[25-27] or of very high-resolution and long inspirals [28, 29]. 
At the same time, the space of parameters is also being care¬ 
fully investigated, both in terms of the variety of the EOSs 
considered [11-13], and of the treatment of radiative losses 
via neutrino cooling [30, 31]. 

When magnetic fields are present, on the other hand, the 


^ In principle, the BMP can also be a stable neutron star, but this would re¬ 
quire that the stars have mass M < /2. Since the mass distribution 

in neutron-star binaries is peaked around 1.3 — 1.4 Mq [24], this implies 
that ^2.8 Mq . Although this cannot be excluded, there is also no 

observational evidence that such massive neutron stars exist. 
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bulk of work carried out so far is considerably more limited. 
Investigations of the impact that magnetic fields have on the 
dynamics of the binary have in fact started with the first short- 
inspiral works in Refs. [32, 33], which were complemented by 
the longer simulations in Ref. [34]. The latter work, together 
with Ref. [35], also investigated the possibility that magnetic 
fields could have an impact on the gravitational-wave signals 
emitted by these systems during the inspiral. The conclusions 
reached were that, for realistic strengths B < 10^^ G, the 
presence of magnetic fields could not be revealed by detectors 
such as advanced LIGO or advanced Virgo. The astrophys- 
ical implications of the merger of a magnetized neutron-star 
binary were explored in Ref. [36], where it was shown that in¬ 
stabilities in the torus orbiting the black hole amplify the mag¬ 
netic field by 3 orders of magnitude and generate a magnetic- 
jet structure characterized by an ordered poloidal magnetic 
field along the black-hole spin axis. The broad consistency 
with the observations in terms of black-hole spin, torus mass 
and accretion rate, and magnetic-field topology offered the 
first evidence that the merger of magnetized neutron stars can 
provide the basic conditions for the central engine of SGRBs. 

Considerable effort has also been dedicated to investigat¬ 
ing the properties of the HMNS under more controlled con¬ 
ditions. For instance, using ultrahigh spatial resolutions but 
axisymmetric initial data. Ref. [37] has provided the first ev¬ 
idence from three-dimensional global simulations that a mag- 
netorotational instability (MRI) [38, 39] is likely to develop 
during the lifetime of the HMNS (see also Refs. [40, 41] for 
earlier work in two dimensions). In addition, again using ax¬ 
isymmetric initial data and different magnetic field configura¬ 
tions, it has been shown that a magnetically driven wind can 
be launched from the outer layers of the HMNS as a result of 
its differential rotation [42, 43]. These works have also high¬ 
lighted that for realistic magnetic field topologies the wind 
is baryon loaded and quasi-isotropic, with bulk velocities of 
^ 0.1 c [43]. More recently, instead, the use of subgrid model¬ 
ing as an effective way to describe the turbulent dynamics that 
develops in the shear layer between the two neutron stars at 
merger, has suggested that amplifications of up to 5 orders of 
magnitude are possible [44], although these amplifications are 
not produced in direct simulations [35, 36], even at very high 
resolution [45]. Finally, progress has taken place also on the 
derivation of improved numerical techniques, such as those 
in Ref. [46], where the significant advantages of a vector- 
potential approach and of a Lorentz gauge were discussed. 

All of the works mentioned above have been carried out 
within the ideal-magnetohydrodynamic (IMHD) approxima¬ 
tion, in which the electrical conductivity is assumed to be in¬ 
finite. Under these conditions, the magnetic fiux is conserved 
and the magnetic field is frozen into the fiuid, being simply 
advected with it. This is a very good approximation for the 
stellar interior before the merger because it neglects any ef¬ 
fect of resistivity on the dynamics of the plasma. After the 
merger, however, there will be spatial regions with very hot 
plasma where the electrical conductivity is finite and the re¬ 
sistive effects, most notably, magnetic reconnection, will take 
place. 

An obvious improvement over the IMHD description is the 


use of the general-relativistic resistive-magnetohydrodynamic 
(RMHD) equations, which provide a complete magnetohy¬ 
drodynamic (MHD) description of regions with a high con¬ 
ductivity, such as the stellar interiors, and of regions with 
small conductivity, such as the electrovacuum exterior. Fur¬ 
thermore, when the conductivity is set to zero, it yields the 
Maxwell equations in vacuum, thus allowing for the study of 
the magnetic-field evolution also well outside the stellar mag¬ 
netosphere [47]. 

Partly because of the increased complexity of the equa¬ 
tions and partly because of the additional difficulties posed 
by their numerical solution (the equations easily become stiff 
in regions of high conductivity), RMHD simulations have 
started only rather recently. Most of the work so far has fo¬ 
cussed on problems in fiat spacetimes [47-53], but general- 
relativistic investigations have also been carried out on fixed 
spacetimes [54]. Indeed, together with the work carried out 
in Refs. [55, 56], those reported here are, to the best of our 
knowledge, the only RMHD simulations of the dynamics of 
binary neutron stars in general relativity. More specifically, 
we have followed the inspiral, merger, and collapse to a black 
hole of a neutron-star binary in which the stars have the same 
gravitational mass of M = 1.625 Mq and are modelled with 
a simple ideal-fluid EOS. 

Complementing the work reported in Refs. [55, 56], which 
concentrated on the electromagnetic emission during the in¬ 
spiral and at the merger, the focus of the simulations reported 
here is that of assessing the impact that resistive effects have 
on the dynamics of the binary after the merger and through 
to collapse to a black hole. To this scope we have carried out 
a close comparison with an equivalent simulation performed 
for the same binary within the traditional IMHD approxima¬ 
tion. In this way it has been possible to determine both the 
similarities between the two regimes and the novel features. 
The most important of such features is the evidence the sur¬ 
vival time of the HMNS before collapse to a black hole in¬ 
creases in a RMHD simulation. This difference is associated 
to a less efficient magnetic braking in the resistive regime, in 
which matter is no longer perfectly advected with the fiow, 
but can move across magnetic-field lines. In turn, this reduces 
the transport of angular momentum away from the central re¬ 
gions of the HMNS, increasing its lifetime. Interestingly, a 
longer-lived magnetized HMNS is of help in those models of 
SGRBs which invoke the existence of a magnetarlike object 
produced after the merger [57-62]. Another important result 
of the simulations reported here, which have been performed 
at higher resolutions and with a different grid structure than 
those in the previous work of Ref. [36], is the confirmation 
that a magnetic-jet structure is formed in the low-density fun¬ 
nel produced by the black hole-torus system. Both in RMHD 
and in IMHD, the magnetic field is predominantly toroidal in 
the highly conducting torus and predominantly poloidal in the 
nearly evacuated funnel. Furthermore, because of the effec¬ 
tive decoupling between the matter and the electromagnetic 
fields achieved in the RMHD simulations, the magnetic-jet 
structure is coherent on the largest scale of our system. How¬ 
ever, as in the IMHD case, also in these RMHD simulations, 
the magnetic-jet structure does not lead to an ultrarelativis- 
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tic outflow. Reconnection processes or neutrino annihilation 
occurring in the funnel could potentially increase the internal 
energy in the funnel and launch a relativistic outflow; none of 
these effects is modeled in our simulations. 

The plan of the paper is as follows. In Sec. II we briefly 
review the mathematical setup of our simulations, concen¬ 
trating mostly on the discussion of the general-relativistic 
RMHD equations used and on the expression of the gener¬ 
alized Ohm’s law we have employed. Section III, on the 
other hand, is dedicated to illustrate the numerical strategy 
employed in the solution of the combined set of the Einstein 
and RMHD equations, including the properties of our compu¬ 
tational grid, of our matching to the low-density exterior, and 
of our initial data. The core of the paper is represented by Sec. 
IV, where we present our results. After a brief overview, we 
discuss the magnetic-fleld topology and magnetic-jet structure 
produced in the simulations, as well as the comparison with 
the IMHD case. Such a comparison goes over a number of as¬ 
pects, from the angular-momentum transfer and lifetime of the 
HMNS, over to the black hole-torus system, and the electro¬ 
magnetic luminosities. Finally, Sec. V contains a conclusive 
summary of our results and the prospects for future research. 
Although this is not the focus of this work, an illustration of 
the magnetic-jet structure obtained in the IMHD simulations 
is presented in Appendix A for completeness. 

We use a spacelike signature (—,+,+,+) and a system of 
units in which c = G = Mq = 1 unless stated differently. 


II. MATHEMATICAL SETUP 
A. General-relativistic RMHD equations 

Much of the numerical setup used in these simulations has 
been presented in greater detail in other papers 118, 34, 63- 
66], and for compactness we will review here only the ba¬ 
sic aspects, referring the interested reader to the papers above 
for additional information. However, given its importance 
here, we will dedicate some space to a review of our fully 
general-relativistic RMHD framework, which was first pre¬ 
sented in Ref. 166] and represents the extension of the special- 
relativistic RMHD formalism discussed in Ref. 147] . A sim¬ 
ilar but independent extension has been presented recently in 
Ref. 160], which describes the first 3 -h 1 general-relativistic 
RMHD implementation in fixed spacetimes. 

We start by presenting the augmented Maxwell equations 

, ( 1 ) 

, ( 2 ) 

where g^^ is the 4-metric, is the Faraday tensor, is 

the Maxwell tensor, is the electric four-current density, and 

0, 7/^ are two auxiliary scalar variables added to the Maxwell 
equations to control the constraints for the magnetic and elec¬ 
tric parts, respectively (see below). 


After a standard 3-^1 splitting of spacetime, the Maxwell 
and Faraday tensors can be decomposed in terms of the elec¬ 
tric {E^) and magnetic fields measured by an observer 
moving along the normal direction (i.e., normal or Eule- 
rian observer) as 

F^^•' = , (3) 

, (4) 

with ;= /y^—g, g the determinant of the 4-metric 

and the Levi-Civita symbol. The same can be done for 

the electric four-current 

:= qn^ + , (5) 

where q and are the charge density and the electric current 
density for an Eulerian observer, respectively. Using these 
definitions and performing a 3-^1 decomposition of Eqs. (1) 
and (2) with respect to the normal vector we arrive at the 
following evolution equations, 

{dt - Cp)E^ - = 

aKE* - aE , (6) 

{dt — + oS/iE'^ = aq — aK, 2 p , (7) 

(dt - Cp)B^ + Vj {aEk) + = aKB ^, (8) 

{dt - CEP + aViB^ = -anp , (9) 

where ^ij is the spatial 3-metric, K := K^- is the trace 
of the extrinsic curvature Kij, a is the lapse and (5 is the 
shift 4-vector. We recall that C (3 denotes the Lie derivative 
along the shift vector and that is related to the four¬ 
dimensional Levi-Civita tensor via or al¬ 
ternatively I s/j, where 7 now is the determinant 

of the 3-metric. The scalar fields 7/; measure the devia¬ 
tion from the constrained solution, with (j) driving the solution 
of Eq. (9) toward the zero-divergence condition = 0, 

and 7p driving the solution of Eq. (7) toward the condition 
ViE'^ = q. This driving is exponentially fast and over a 
time scale 1/k,. This approach, named hyperbolic divergence 
cleaning in the context of IMHD, was proposed in Ref. 167] 
as a simple way of solving the Maxwell equations and enforc¬ 
ing the conservation of the divergence-free condition for the 
magnetic held. This method has been extended to the resistive 
relativistic case in Refs. 147, 49]. 

An obvious consequence of the Maxwell equations in 
RMHD is the conservation law associated with the electric 
charge 

V^/^ = 0, (10) 

which provides an evolution equation for the charge density 

{dt-C^)qE\/i{aE)=aKq. (11) 

Combining the MHD and Maxwell equations we obtain the 
following set of evolution equations, which we write in a flux- 
conservative form as 
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, (12a) 

dt{^E^) + dk{-p^^E^ - ae^’^^y/jBj) = -^E\dkP^) - - a^E , (12b) 

dt4> + 4(-/3V + aS'') = -4>{dkP'^) + B>^{dka) - - aK<j>, (12c) 

dti^ + 4(-/3V + ctE’^) = -i’idkP'^) + E'^idka) - |(7""47im)^'' + aq - an^P, (12d) 

dt{Vjq) + dk[^/j{-l3^q + aj'^)]=0, (12e) 

dt{ViD) + dk[Vi{-l3’‘D + av’^D)]=0, (12f) 

dtiVir) + dkiVli-fi’^T + aiS’^ - v'^D)]} = ^{aS^^Ki^ - S’^d^a), (12g) 

dtiViSi) + dkiVli-l^’^Si + aS’^^] = V7 [^S^^'daim + S^di/S’^ - (r + D)5,a] . (12h) 


The fluid variables U, Si, and Sij are the conserved rest-mass density, the conserved energy density, the conserved momen¬ 
tum, and the fully spatial projection of the energy-momentum tensor, respectively. Their explicit deflnitions are therefore 


D:=pW, (13) 

T:=U-D = phW‘^-p+^{E‘^ + B‘^)-D, (14) 

Si := phW\ + CijkE^B ’^, (15) 

Sij := phW‘^ViVj -h ^ijP - EiEj - BiBj + ^jij{E‘^ + B‘^) . (16) 


Here, IT^ = au^ = Ui/vi is the Lorentz factor, where are 
the components of the fluid 4-velocity and v'^ are the compo¬ 
nents of the 3-velocity as measured by the Eulerian observer. 
Furthermore, h := (e+p )/p = 1+e+p /p is the enthalpy, with 
e = p(l + e) the total energy density, p the pressure, e the spe- 
ciflc internal energy, and p the rest-mass density [68]. The im¬ 
portant difference between the RMHD and IMHD equations 
is that they involve stiff relaxation terms that pose serious nu¬ 
merical limitations on the time evolution of the equations. For 
this reason, a distinct class of implicit-explicit evolution meth¬ 
ods has been developed, the RKIMEX schemes [69], which 
we have presented in detail in Ref. [66]. 

1. Generalized Ohm’s law 

To close the system of equations presented above, a rela¬ 
tion for the electric current density in terms of the other flelds 
is necessary, just like Ohm’s law provides a prescription for 
the spatial conduction current to be proportional to the elec¬ 
tric held. A generalized Ohm’s law provides the necessary 
coupling of the electric current density to the electromagnetic 
and matter flelds. Previous work toward relativistic versions 
of the generalized Ohm’s law includes the investigations in 
Refs. [70-73]. In one of the simpler cases the spatial conduc¬ 
tion current can be considered as proportional to the electric 
held measured by the comoving observer. Therefore the elec¬ 
tric four-current density can be written as the superposition of 
an advective and a conductive current [74], which takes the 
form of the generalized Ohm’s law 

= + = + ( 17 ) 


Here, q := —I^u^ = + (Jan^)]/IF is the electric charge 

density measured in the rest frame comoving with the fluid, 
and should be contrasted with q := which is instead 

the charge density measured by the Eulerian observer. Sim¬ 
ilarly, and are, respectively, the electric held and the 
electrical conductivity of the medium (which is a rank-2 sym¬ 
metric tensor) as measured in the same frame. In collisional 
plasmas the current in the comoving frame can be considered 
to be carried by the mobile electrons, with charge e, and the 
conductivity tensor becomes 

, (18) 

with 6^ being the magnetic held in the comoving frame, ^ := 
CTe/me, (T := nee^/(l+6^), andre the electron collision 
time scale. 

Expressing the four-current density of Eq. (17) in terms of 
the flelds measured by the Eulerian observer we arrive at 

+ Wa[E^ + - {v^E^)v^]^ 

Waf{E°‘Ba)[B^^ - (19) 

In deriving Eq. (17) we have made the implicit assumption 
that the collision frequency between particles is much larger 
than the typical oscillation frequency of the plasma, which, we 
recall, is deflned as cCp := (dyrnge^/me)^/^. This implies that 
electrons and ions can reach equilibrium on very short time 
scales and any correction due to the mass difference between 
electrons and protons can be neglected. As a result, there is no 
global charge separation and the plasma is neutral. Note also 
that the first term in Eq. (18) accounts for an isotropic scalar 
law for the current, while the rest represent anisotropies due 
to the presence of a magnetic held in the comoving frame. 
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Ideally, it would be desirable to have a well-defined pre¬ 
scription of the conductivity tensor as a function of the fiuid 
properties, which stems from the mi¬ 

crophysical properties of the plasma (see Ref. [73] for a re¬ 
cent discussion). In practice, however, we are far from having 
such a prescription in the extreme physical conditions char¬ 
acterizing merging neutron stars. However, if the collision 
time scale is much smaller than the electron cyclotron pe¬ 
riod^ or, equivalently, if ^ ^ electrons do not have 

sufficient time to gyrate perpendicular to the magnetic field 
lines. Under these conditions, the isotropic part of the con¬ 
ductivity is the dominant one (electrons essentially slide along 
the magnetic-field lines), and expression (18) can be approx¬ 
imated as ^ cfQhv As a result, the spatial three-current 
density in the Eulerian frame coming from the generalized 
Ohm’s law (19) can be simplified so that 

J* = qv^ + Wa[W + - {vkE’^Y ], (20) 

where we used the fact that the four-current density can be 
also written as /^ = qn^ + J^. In our simulations, the con¬ 
ductivity (7 is chosen to be either a constant or a function of 
the rest-mass density and a discussion will be presented in 
Sec. me 1. The last term in Eq. (19) represents the Hall ef¬ 
fects, which, however, we set to zero for simplicity. 


III. NUMERICAL SETUP 

A. Eield equations 

The evolution of the spacetime (i.e., of the 3-metric, ex¬ 
trinsic curvature, and conformal factor) is obtained using the 
Me La chi an code, which implements the BSSNOK formu¬ 
lation of the Einstein equations [75-77] employing three- 
dimensional finite-differencing operators for calculating the 
fiuxes on the right-hand side of the Einstein equations [78]. 
The time integration is carried out using the third-order ac¬ 
curate strong stability preserving implicit-explicit Runge- 
Kutta scheme outlined in Ref. [66]. The time step on each 
grid is limited by the Courant-Eriedrichs-Lewy (CEL) condi¬ 
tion [68] and hence the Courant coefficient is set to be 0.25 
on all refinement levels. A Kreiss-Oliger type dissipation is 
added to the spacetime evolution equations to ensure that any 
high-frequency noise produced during the evolution (mostly 
at the refinement levels boundaries) is damped. The damping 
parameter in the evolution equation for the shift was carefully 
chosen to be 0.71 so that the additional term does not become 
stiff on the coarser grids. 


^ The electron cyclotron period is defined as Pc,e ’•= 27r/a;c,e, where 
<^c,e := eBlrrie is the cyclotron frequency and represents the frequency 
at which electrons gyrate perpendicular to the magnetic-field lines. 


B. GR-RMHD equations 

The general-relativistic RMHD (GR-RMHD) equations are 
solved with the WhiskyRMHD code [66], which uses high- 
resolution shock-capturing schemes [68, 79]. In particular, the 
reconstruction of the conserved variables is achieved via the 
Piecewise Parabolic Method [80], while the fiuxes are calcu¬ 
lated through the approximate Riemann solver introduced by 
Harten-Lax-van Leer-Einfeldt [81] and which requires only 
knowledge about the maximum characteristic speeds of the 
system, i.e., the speed of light in this case. 

The system of RMHD equations is closed by describing the 
fiuid as ideal and with a T-law EOS 

p = pe{T -l), (21) 

where T = 2 and K = 123.6. This EOS is clearly an ide¬ 
alization. Much more sophisticated EOSs have been imple¬ 
mented in our numerical infrastructure [13, 82, 83] and have 
been used by other groups in general-relativistic magnetohy¬ 
drodynamic (GRMHD) simulations [31]. However, this ideal¬ 
ization is probably adequate at this stage as here we are mostly 
interested in assessing the differences introduced in the dy¬ 
namics of the BMP by resistive effects. These effects suffer 
from even larger uncertainties than those associated to the dif¬ 
ferent EOSs. 


C. Adaptive mesh refinement and symmetries 

Our code makes use of the Cactus [84] computational 
framework, which allows us to employ a box-in-box vertex- 
centered adaptive mesh refinement grid hierarchy that tracks 
the “center of mass” of the stars as they orbit each other. This 
was achieved via the BNSTracker thorn implemented by W. 
Kastaun and the Carpet driver [85]. The numerical domain 
consists of six levels of refinement with the resolution dou¬ 
bling between adjacent refinements. In addition, we employ 
moving refinement boxes in order to track the high-density 
regions. The outer boundary of the computational domain is 
located at « 378 km. The finest resolution during the inspi¬ 
ral and merger is Ax « 296 m, but an extra refinement level 
is activated right after the merger with a resolution of Ax 
148 m. It is important to remark that, although in a small re¬ 
gion of spatial extent ^ 13.5 km, this resolution is consider¬ 
ably higher than the one used in Ref. [36] (i.e., « 221 m with 
a spatial extent of ^ 35.4 km), where the first evidence was 
given that the merger of a binary system of magnetized neu¬ 
tron stars can lead to the formation of a magnetic-jet structure. 
Hence, although in a resistive framework, the simulations re¬ 
ported here can be considered as a “higher-resolution” coun¬ 
terpart of those presented in Ref. [36]. We also note that when 
considering initial data consisting of a binary in quasicircu¬ 
lar orbits (see Sec. IV C 2) we do not activate the extra re¬ 
finement level after the merger, hence keeping a resolution of 
Ax 296 m to describe the HMNS. We do this so as to have 
two different resolutions to describe the HMNS and hence to 
study the numerical consistency of the solution. 
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To reduce the computational cost associated with the nu¬ 
merical evolution of an equal-mass binary, we use a reflection- 
symmetry condition across the ^ = 0 plane and a 7r-symmetry 
condition across the x = 0 plane. These conditions respect 
the symmetries of the scenario we are investigating. The outer 
boundary conditions are set by using simple zeroth-order ex¬ 
trapolation of the hydrodynamic variables. For the electro¬ 
magnetic and spacetime variables, we employ simple Som- 
merfeld radiative boundary conditions because of the nature 
of the flelds. Given the long time scale over which our simu¬ 
lations are carried out, reflections due to imperfections of the 
boundary conditions experience several domain crossings be¬ 
fore the end of the simulation. We plan to improve on this in 
the future by the application of maximally dissipative bound¬ 
ary conditions that minimize the effects of reflections. 


7. Exterior matching and atmosphere handling 


Our goal with the use of a RMHD framework is that of 
modelling the exterior of the neutron stars as an electrovac¬ 
uum, where both the conductivity and the charge density are 
negligibly small, so that electromagnetic flelds should obey 
in these regions the Maxwell equations in vacuum. On the 
other hand, we want to model the interior of the stars as highly 
conducting, so that our equations recover the IMHD limit in 
such regions. There are several different ways to achieve this; 
see, for example. Refs. [66, 86-88]. Each of them has, in 
our view, its advantages and disadvantages. However, because 
they all try to model the difficult transition region between two 
regimes that are intrinsically different, they all represent an 
approximation. This difficulty is not speciflc to this problem 
but is a typical feature of physical problems as, for example, 
in the transition from an optically thick to an optically thin 
regime in radiative-transfer calculations. Our approach is also 
an approximation and is similar to the one in Refs. [66, 87] 
in the sense that the matching from the stellar interior to the 
stellar exterior is achieved through a carefully chosen con¬ 
ductivity profile. More specifically, the conductivity profile 
adopted is directly related to the conserved rest-mass density 
(and hence to the rest-mass density) and given by 


cr := (Jo max 


1 + exp [2Dto\{D - 

-^rel)/ -^atm] 


0 , 
( 22 ) 


where do ^ 1 corresponds to a uniform scalar conductivity. 
The parameters T)toh and determine how sharp the tran¬ 
sition to the exterior is. For the simulations reported here, we 
have chosen^, ctq = 10^ = 2.0 x 10^^ sec“^, T)toi = 0-01 
and D^eX — 100 -Datm — 100 Patm? where -Datm and Patm are 


the values of the conserved and primitive rest-mass density in 
the “atmosphere.” 

We recall that, as in other Eulerian hydrodynamics and 
MHD codes, also our WhiskyRMHD code makes use of a 
very low rest-mass density fluid to handle the evolution of the 
MHD equations in regions which are associated to the exte¬ 
rior of the stars. In such a region, we follow the same pre¬ 
scription initially implemented in Ref. [89] and then adopted 
in essentially all of the simulations performed with our code, 
even in its IMHD incarnation [63]. In essence, we treat as 
atmosphere any region in the computational domain which is 
below Patm, which we take here to be 6.17 x 10^ g cm“^, that 
is approximately 8 orders of magnitude smaller than the max¬ 
imum rest-mass density."^ In such a region, we set the fluid 
3-velocity to zero, the rest-mass density to a floor value, and 
the speciflc internal energy to the value it assumes for a fluid 
following a polytropic EOS [68].^ In addition, and differently 
from the IMHD implementation of Ref. [35], the use of the 
prescription (22) automatically sets the conductivity in the at¬ 
mosphere to zero, so that the Maxwell equations reduce there 
to the Maxwell equations in vacuum. Hence, in our prescrip¬ 
tion the atmosphere is de facto a cold, static, uniform fluid in 
which electromagnetic waves propagate as if in vacuum. 

A few remarks should be made about modelling the neu¬ 
tron stars’ exterior, which in reality is expected to be a 
highly conducting, low-density, possibly magnetically dom¬ 
inated plasma in very strong magnetic flelds. 

First, with our prescription for the atmosphere and with the 
conductivity (22), the latter is zero also in regions that are 
not at the atmosphere level but close to it, i.e., at D < Drei- 
We do this because the strong winds that are produced at the 
merger and later on rapidly All the computational domain. A 
direct consequence of this is that the atmosphere is present 
only close to the boundaries and therefore an IMHD prescrip¬ 
tion would be met essentially everywhere. Yet, we are inter¬ 
ested in deviations from perfect-flux freezing and in particular 
in regions where plasma is tenuous and the magnetic flelds es¬ 
sentially decouple from the matter. We can effectively achieve 
this by setting cr = 0 in any region that is below Drei- This ap¬ 
proach therefore accomplishes our goal of decoupling in the 
low-density regions the evolution of the electromagnetic flelds 
from the dynamics of the plasma. Of course, an electrovac¬ 
uum prescription where the rest-mass density is nonzero is, 
strictly speaking, inconsistent, but we believe that this is a tol¬ 
erable inconsistency, given that the rest-mass density in these 
regions takes essentially the smallest values in the whole do¬ 
main. 

Second, albeit somewhat arbitrary, our approach suffers 
from the same uncertainties of other approaches suggested 
in the literature to match the two different regimes (see, 
e.g.. Refs. [56, 87, 90]). Ideally, its robustness can be vali¬ 
dated by varying the free parameters Drei and Patm, although 


^ We recall that in units in which the speed of light is not set to 1, : = 

L^(7o/c^ represents the Ohmic diffusion time scale, with L the typical 
length scale of the field. For L ~ 10® cm, ~ 10^ sec, which is obvi¬ 
ously much larger than the timescale over which our simulations are carried 
out, i.e., ~ 10“^ sec. 


^ In practice, to avoid being sensitive to the threshold value, we set to atmo¬ 
sphere any cell of which the rest-mass density is below patm + Ptoh where 
Ptol ~ 10 ^ Patm [18]. 

^ A different prescription is used for the specific internal energy in the case 
of hot, nuclear-physics EOSs (see Ref. [82] for details). 
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this is something that admittedly we are not able to do here be¬ 
cause of the computational costs involved. More importantly, 
we find that this approach allows us to take an important step 
beyond the previous IMHD treatment presented in Ref. [91]. 

Third, while matching to a force-free regime is appealing, 
it can be rather dangerous in the physical conditions encoun¬ 
tered after the merger (by contrast, the force-free approxi¬ 
mation is probably very good before the merger or in black 
hole-neutron-star binaries, when only magneto spheric effects 
are expected to take place). Our calculations reveal in fact 
that the exterior of the conducting matter (either the HMNS 
or the torus) is always matter dominated and the plasma beta 
parameter, i.e., the ratio of the gas-to-magnetic pressure,^ 
:= 2plB‘^, is at least 10^ and of the order of 10^ in the 
polar regions. This is far from the condition of 1, 

where the matter inertia can be neglected and the force-free 
approximation is a good one. 

Finally, despite the fact that our implementation offers a 
control over the amount of resistivity in different parts of the 
fiow, the choice of realistic values for the resistivity is far from 
trivial (see, e.g.. Ref. [92]) and not addressed at all in these 
simulations. 


2. Miscellanea 

As already mentioned in Sec. II A, in order to ensure that 
the magnetic field is essentially divergence free, we have em¬ 
ployed the divergence cleaning method of Ref. [67], which, 
however, requires choosing a suitable value for the constant 
n [cf., Eqs. (7) and (9)]. In the simulations presented here, 
we have set n = 0.075. This value does not lead to stiffness 
problems in the coarser grids and at the same time provides 
a rapid damping of the constraint violation on a time scale 

When after the merger the HMNS collapses to a black hole, 
steep gradients appear in the rest-mass density and cover just 
few grid points of the finest grid. If the resolution is suffi¬ 
ciently low, these gradients are simply dissipated numerically 
and the evolution can proceed without problems. However, 
for the rather high resolutions used here for the finest grid, 
i.e., 148 m, the numerical dissipation is smaller and the gra¬ 

dients are not removed, leading to failures in the conversion 
of the conserved variables to the primitive ones. To counter 
this problem we reset to atmosphere the hydrodynamical vari¬ 
ables in a mask inside which the lapse function goes below a 
threshold, e.g., a < athr = 0.1. This reset is done only for 
the hydrodynamical variables, while the spacetime and elec¬ 
tromagnetic ones are evolved as usual. Furthermore, the reset 
is applied in practice to a handful of cells, well inside the ap¬ 
parent horizon and thus not influencing the matter dynamics. 


^ We note that our definition of /3p is the one normally used in plasma 
physics, but the inverse of the one employed in other numerical-relativity 
calculations, e.g.. Ref. [21]. 


D. Initial data 

The initial data consist of a magnetized binary neutron- 
star system of total Amowitt-Deser-Misner (ADM) mass 
M^dm =3.25 Mq and an initial orbital separation of 45 km. 
Each star has a baryon mass equal to 1.625 Mq and an equa¬ 
torial radius of i?eq = 13.68 km, so that the initial separation 
corresponds to approximately 3.3 i?eq. The initial orbital ve¬ 
locity is flo = 1-85 rad ms“^, and the maximum rest-mass 
density is 5.91 x lO^^grcm”^. Lacking self-consistent ini¬ 
tial data for magnetized binaries, our initial data are gener¬ 
ated by the LORENE library as an unmagnetized irrotational 
binary in equilibrium on a quasicircular orbit [93, 94]. The 
magnetic field is then superimposed on the unmagnetized 
constraint-satisfying solution. Following Refs. [34, 35], the 
initial magnetic field is fully contained inside the stars and 
purely poloidal. This is achieved after prescribing the toroidal 
vector potential to have the form 

= Ab [inax(P - Pcm, 0)]^ , (23) 

where Pcut = 0.04 Pmax determines the point at which the 
magnetic field goes to zero (typically before it reaches the 
surface). The resulting (poloidal) magnetic field is just the 
curl of the vector potential and leads to a maximum magnetic- 
field strength of 1.97 x 10^^ G at the pressure maximum. Be¬ 
cause the initial data constructed in this way are not a solution 
of the full Einstein-Euler-Maxwell system, they will intro¬ 
duce an increased violation of the constraint equations, which 
amount to ^ 2.5 x and ^ 1.7 x in 

the P 2 -norm, of the Hamiltonian and momentum constraints, 
respectively. These values should be compared with those 
before the introduction of the magnetic field and which are 
about 1 order of magnitude smaller, i.e., 5 x 10 
1.5 X 

The perturbations introduced by the addition of the mag¬ 
netic fields are small enough so as not to have a significant 
effect on the dynamics of the binary. Indeed, our experience 
is that the P 2 -norm of the constraints relaxes to values com¬ 
parable to those of simulations of unmagnetized binaries after 
about one crossing time or, equivalently, one orbit. Obvious 
ways to improve this approach, and which could become im¬ 
portant for an accurate modelling of the gravitational-wave 
emission during the inspiral, exist. Among them are the use 
of consistent initial data for magnetized binaries or the sim¬ 
ulation of binaries with much larger separations, so that the 
system has several orbits to reach a more consistent MHD 
equilibrium. 

In addition, because the main focus of this work is the as¬ 
sessment of resistive effects on the postmerger dynamics and 
hence a comparison between IMHD and RMHD simulations, 
our interest in an accurate treatment of the initial data is rather 
limited here. As a result, and once again to reduce the compu¬ 
tational costs, we accelerate the inspiral as first suggested in 
Ref. [95]. More specifically, for most of our simulations we 
modify the initial linear momenta by adding an initial inward 
radial velocity which is ^ 20% of the orbital velocity. This 
reduces the number of orbits at this separation from ^ 3.5 to 
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only ^ 1.5. Of course, the constraint violations introduced 
in this way are even larger than those discussed above with 
the introduction of magnetic fields^, but we have also verified 
that this does not introduce qualitative differences by compar¬ 
ing the results of these grazing collisions with those obtained 
from the corresponding binaries in quasicircular orbits. 

IV. RESULTS 

The dynamics of the same neutron-star binary considered 
here has been previously investigated in hydrodynamic sim¬ 
ulations [18], in IMHD simulations [35], and, more recently, 
also in RMHD simulations [56]. In all cases, it was shown 
that after the merger this specific model forms a rapidly ro¬ 
tating HMNS with a high degree of differential rotation. The 
transient object collapses later to form a rotating black hole of 
mass ^ 2.9 — 3.0 Mq. The magnetic-field strengths chosen 
here and in previous works [35, 56] are not sufficiently high to 
affect the bulk dynamics of the binary during the inspiral [34], 
which can be considered to be equivalent to the purely hydro- 
dynamical case for all practical purposes. 

In the following subsections, we focus on the results ob¬ 
tained from the application of our RMHD implementation. 
We start by providing a general description of the basic fea¬ 
tures of the RMHD dynamics (Sec. IV A) and then move to 
make a detailed comparison with an IMHD simulation of the 
same neutron-star binary (Sec. IV C). We note that our focus 
here is different from the one in Ref. [56], which matched the 
resistive description to a force-free one to study the interaction 
of the two stellar magnetospheres before the merger. Here, on 
the other hand, we are mostly interested in the postmerger ob¬ 
ject and on the effects that resistivity has on the dynamics of 
the HMNS and subsequent black hole-torus system. Because 
of this, and because the scenario we are investigating is pol¬ 
luted by the large baryonic winds produced after merger, our 
resistive matching is made to an electrovacuum exterior. In 
this sense, the work carried out here and in Ref. [56] provide 
a complementary description of the dynamics of magnetized 
binary neutron stars in full general relativity. 

A. Rapid overview 

The dynamics of the binary when evolved within the 
RMHD framework is summarized in Figs. 1 and 3. The first 
one, in particular, reports two-dimensional cuts in the {x^y) 
and (x, z) planes of the rest-mass density p (left panel), of the 
specific internal energy e (middle panel), and of the modulus 
of the magnetic field \B\ := (5*5^)^/^ (right panel). Marked 


^ For completeness, we can compare the violations in the L 2 -norm of the 
Hamiltonian and momentum constraints for binaries with reduced initial 
momenta with the corresponding violations for binaries on quasicircular 
orbits. The latter amount to ~ 2 x 10“for the Hamiltonian 

/ ADM 

constraint and to ~ 4 x average of the momentum 

constraints. 


with white lines are the projection of the magnetic-field lines 
on the different planes, while marked with green solid lines 
are the isocontours of the rest-mass density. The snapshots 
refer to times t = 2.43 ms (top row), t = 3.91 ms (middle 
row), and t = 11.60 ms (bottom row). The positions of the 
two stars are marked with x and + symbols, and the stars’ 
trajectories are marked with red dashed lines. 

Given the small initial orbital separation of 45 km and the 
reduced linear moment, the two stars merge very rapidly. 
More specifically, at approximately t 0.5 ms, the two 
stellar surfaces start entering in contact, although the actual 
merger takes place at t ^ 3.91 ms.^ As the merger takes 
place, the two stellar cores become significantly distorted by 
the large tidal fields and produce spiral arms. At the lead¬ 
ing edges of these spiral arms, the specific internal energy in¬ 
creases through shock heating (cf., the middle column of Fig 
1 ). 

The time t = 2.42 ms in the top row of Fig. 1 corresponds 
to one orbit of the binary, which is sufficient for the magnetic 
field to diffuse over the thin transition layer close to the sur¬ 
face of the stars in an attempt to settle to a new equilibrium 
configuration (cf., the discussion in Sec. IV.C.l of Ref. [66]). 
Once the magnetic field has diffused out of the star it contin¬ 
ues to propagate also in regions that are dynamically treated 
as atmosphere and where the electrical conductivity is set to 
zero as if the medium was an electrovacuum. In this way, 
we achieve a rather smooth transition between the highly con¬ 
ducting stellar interior and the electrovacuum exterior. This 
is shown in Fig. 2, where we show the electrical conductivity 
at three reference times on the (x, z) plane. Note that, the re¬ 
gion in black corresponds to our electrovacuum but does not 
coincide with the atmosphere. Indeed, the region in black is 
filled with tenuous plasma as can be seen in the left columns 
of Fig. 1. Note also that at this time the magnetic-field topol¬ 
ogy is still predominantly poloidal in the exterior of the star. 
However, a toroidal component is also being generated as the 
highly conducting material in the stellar interior shears the 
poloidal magnetic field in the lower-density, high-conductivity 
spiral arms. 

As mentioned above, at t = 3.91 ms the merger takes 
place, at least as measured from the position of the first peak 
in the gravitational-wave amplitude. When this happens, a 
vortex sheet is created between the two stars, which could 
lead to the onset of a Kelvin-Helmholtz instability and the 
generation of a large-scale and ultrastrong magnetic field [97] 
(see the middle row of Fig. 1). It is presently a matter of de¬ 
bate whether such large-scale magnetic fields can be produced 
with amplifications of several orders of magnitude. Present 
direct simulations are not able to reach the resolutions nec¬ 
essary to resolve the turbulent motion produced by the insta¬ 
bility [34, 35]. The results obtained so far with direct, very 
high-resolution simulations, either local [98] or global [45], 
indicate that the amplification of the magnetic field is of a fac¬ 
tor 20 at most, most likely because resistive instabilities dis- 


^ As is customary, we define the time of merger as the time of the first peak 
of the gravitational-wave amplitude [13, 83, 96]. 
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FIG. 1. Snapshots of two-dimensional cuts in the (x,y) and {x,z) planes of the rest-mass density p (left panel), of the specific internal 
energy e (middle panel), and of the modulus of the magnetic field \B\ := (right panel). From the top, the snapshots correspond 

to times t = 2.43 ms (top row), t = 3.91 ms (middle row), and t = 11.60 ms (bottom row). Shown with white lines are the projection 
of the magnetic-field lines on the different planes, while marked with green solid lines are the isocontours of the rest-mass density at p = 
{6.2x10'^, 1.2xl0^\ 2.5xl0^\ 3.7xl0^\ 4.9xl0^\ 6.2xl0“, 1.3xl0^^ 2.5xl0^^ 3.7xl0^^ 4.9xl0^^ 6.2 x 10^^} gr cm“^ 
The positions of the two stars are marked with x and + symbols, and the stars’ trajectories are marked with red dashed lines. The different 
panels represent the evolution of the HMNS and highlight that no ordered magnetic-field topology emerges; this will change when the HMNS 
collapses to a black hole (see also Figs. 3-6). 
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FIG. 2. Two-dimensional cuts on the (x^ z) plane of the electrical conductivity distribution at different times t — 
{11.601, 16.141, 18.537} ms (left to right). Note that, the conductivity is very large in the HMNS (left panel) and torus (middle and 
right panel), where the IMHD is recovered. On the other hand, the conductivity is outside the HMNS/torus, where the electrovacuum limit is 
reached. Note also that the region in black does not coincide with the atmosphere, but is filled with tenuous as can be seen in the left columns 
of Figs. 1 and 3. 



FIG. 3. The same as in Fig. 1 but for times t = 15.61 ms (top row) and t = 18.54 ms (bottom row), when a black hole has already been 
formed. Note that, in contrast to the dynamics of the HMNS shown in Fig. 1, the collapse of the HMNS leads to the generation of large-scale 
coherent magnetic fields and the emergence of a magnetic-jet structure around the black-hole rotation axis. The magnetic field is mostly 
poloidal in the funnel and mostly toroidal in the torus (see also Figs. 5 and 6). 
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FIG. 4. Two-dimensional cuts on the (x, z) plane of the specific internal energy and of the velocity field in this plane. Note that the times 
selected coincide with those presented in Fig. 2, and are useful to gain information on the properties of the flow during the HMNS stage and 
after the collapse to a black hole. The flow is essentially given by the magnetically driven wind during the HMNS stage (left panel). However, 
once a funnel is produced, the flow is ingoing at heights < 100 km, but becomes outgoing for z > 100 km (middle and right panels). 


rupt the Kelvin-Helmholtz unstable vortex [98] (but see also 
Ref. [44] for recent simulations with subgrid modeling which 
could lead to much larger amplifications). 

In the simulations reported here, we find that the magnetic- 
field magnitude increases slightly less than an order of mag¬ 
nitude during this stage (see Fig. 12 and the discussion in 
Sec. IV C). This moderate growth of the magnetic field is pos¬ 
sibly due to insufficient resolution and our inability to capture 
the dynamics of the relevant scales. On the other hand, it could 
also indicate that a Kelvin-Helmholtz instability simply does 
not develop. A possible reason is that the time scale of the in¬ 
stability might be longer than the dynamical time scale and the 
shearing motion could very rapidly get destroyed as the two 
stellar cores collide on a time scale of a fraction of a millisec¬ 
ond. We should also note that resistive effects might become 
important at this stage as the reconnection of magnetic-field 
lines might lead to the acceleration of matter due to Ohmic 
heating. We believe that higher-order schemes, as those pre¬ 
sented in Refs. [25, 26, 48], could help resolve this issue. 

The bar-deformed HMNS produced after the merger is dif¬ 
ferentially rotating, and magnetic braking transfers angular 
momentum from the inner core to the outer parts of the star. 
The spiral arms widen and merge together generating more 
shock heating and dissipation. A magnetically driven wind 
as a result of differential rotation is launched from the outer 
layers of the HMNS [42, 43]. The wind could play an impor¬ 
tant role in the modelling of short gamma-ray bursts (SGRBs), 
which show an extended x-ray emission [61, 62]. The wind 
is not constant in time, but rather characterized by a bursty 
activity in which high internal energy plasma blobs (i.e., lo¬ 
cal concentrations of specific internal energy a few kilometers 
in size) are launched from near the black-hole horizon and 
propagate along the ^-direction (bottom row of Fig. 1). In¬ 
terestingly, the bursts observed in the specific internal energy 
of the hot rotating halo that forms around the central object 
are anticorrelated with the bursts observed in the modulus of 
the magnetic field; a more detailed discussion on these bursts 
follows in Sec. IV C 4. 

Unfortunately, the resolution used here is insufficient to be 
able to track the development of an MRI, which is, however, 
expected to develop [37, 45] and could significantly amplify 
the magnetic field. Blind to this effect, our simulation shows 
that at t = 11.60 ms (see the bottom row of Fig. 3) magnetic 
braking has managed to store enough rotational energy in the 


winding of the magnetic-field lines so that the inner core of the 
star is now less differentially rotating. The direct consequence 
of this is that the HMNS collapses to a black hole of mass 
M = 2.88 Mq and dimensionless spin a = J/M‘^ = 0.87, 
as measured from the apparent horizon [99, 100]. We post¬ 
pone the discussion on how angular momentum is transported 
outward and how this affects the lifetime of the HMNS to 
Sec. IV C 2. 

At time t = 15.61 ms, the black hole is surrounded by 
a thick accretion torus that is responsible for confining and 
collimating the magnetic field along the z-axis (cf., the top 
row of Fig. 3). The properties of the black hole-torus system 
are shown in Table I at approximately 4.74 ms after the col¬ 
lapse and show that the mass of the torus is O.OQSMq. The 
magnetic-field topology and matter dynamics soon after the 
collapse are highly turbulent, but the high degree of symme¬ 
try introduced by the black hole, which is gravitationally dom¬ 
inant over the torus, rapidly establishes some order in this sys¬ 
tem. After about one orbital period, in fact, the torus becomes 
essentially axisymmetric and the matter in the polar region is 
rapidly accreted onto the black hole, giving rise to a funnel 
where the rest-mass density reaches values close to that of the 
atmosphere (cf., the bottom row of Fig. 3). 

After black-hole formation, the plasma dynamics in the 
funnel is far from being stationary and continues with repeated 
bursts having a period of about 2.4—3.7 ms. While a more de¬ 
tailed discussion of these bursts is postponed to Sec. IV C 4, it 
is useful to remark here that the ejected material does not have 
sufficient energy to reach large distances away from the black 
hole. This is probably due to the fact that, although the mag¬ 
netic field is comparatively strong, the material in the funnel is 
still matter dominated, with r\j 10^ — 10^. Such large val¬ 
ues are not particularly surprising since the initial magnetic 
field in the stars is rather small, i.e., ^ 10^^ G, and is not 
amplified significantly. At the same time, the torus angular 
velocity profiles have become nearly Keplerian, and the MRI 
could develop (cf.. Fig. 14). However, as for the HMNS, also 
the spatial resolution of the grid covering the torus is too small 
to capture the fastest growing modes of the instability in the 
torus (see discussion in Sec. IV C 5). 

For completeness, we report in Fig. 4 three different two- 
dimensional cuts on the (x, z) plane of the specific internal 
energy and of the velocity field in this plane. The times se¬ 
lected are the same as those presented in Fig. 2, and are use- 
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FIG. 5. Large-scale two-dimensional snapshots on the {x, z) planes of the rest-mass density (top row) and of the magnetic field (bottom row). 
The two columns refer to t = 11.60 ms (left column), when the HMNS has not yet collapsed and to t = 18.54 ms (right column), when a 
black hole has already been formed. Note again the formation of a magnetic-jet structure around the black-hole rotation axis, which extends 
on scales that are much larger than those of the accreting torus (i.e., ~ ±60 km) and of the black hole (i.e., ~ ±5 km) (see also Figs. 3, 6, 16, 
and 17). 


ful to gain information on the properties of the flow during 
the HMNS stage and after the collapse to a black hole. Note 
that during the HMNS stage the flow consists mostly of the 
intense magnetically driven outgoing wind discussed above 
(left panel). However, once a funnel is produced, the flow 
is ingoing at heights z < 100 km but becomes outgoing for 
^ > 100 km (middle and right panels). The location of this 
stagnation point varies with time and follows the expansion 
of hot fluid which is produced by the accretion process and 
that can be followed via the increases in the speciflc internal 
energy. In all cases considered, the vertical flow very close to 
the black hole is ingoing. 


B. Magnetic-field topology and magnetic-jet structure 

As mentioned above, the high degree of symmetry near 
the rapidly rotating black hole induces a quick rearrangement 
of the matter and of the magnetic fields. As a result, the 
magnetic-field topology at time t = 18.54 ms changes in the 
funnel and develops a dominant poloidal component, giving 
rise to a well-defined magnetic-jet structure, which is almost 
axisymmetric. This result is similar to what was already found 
in the simulations of Ref. [36], with the important difference 
that this configuration has been reached with a higher spatial 
resolution and a consistent treatment of the resistivity. 

The rest-mass density of the plasma in the funnel is close to 
that of the atmosphere, but also slightly larger. Hence, given 
our choice of the conductivity profile in Eq. (22), the conduc¬ 


tivity is essentially zero everywhere in the funnel (see the mid¬ 
dle and right panels of Fig. 2). This has two important conse¬ 
quences. First, the dynamics of the electromagnetic fields in 
this region is not that prescribed by the IMHD equations but 
rather that of electromagnetic waves in vacuum. At the same 
time, because of its (comparatively) small rest-mass density 
and pressure, the matter in the funnel tends to move along 
the field lines. We should clarify that this behavior is not 
achieved because the test-particle limit of the RMHD equa¬ 
tions is reached (the rest-mass density and pressure are in fact 
nonzero) but rather because the matter in the funnel can only 
move in the vertical direction, either accreting onto the black 
hole or moving outward (the matter in the funnel has low or 
zero specific angular momentum). 

Modelling the dynamics of the matter in the funnel is 
among the most challenging aspects of these calculations. Al¬ 
though the matter there has the largest magnetic-pressure sup¬ 
port (i.e., > 10^), this is still about 4 orders of magnitude 

away from being magnetically dominated. The reason for this 
behavior is most likely due to the comparatively weak mag¬ 
netic fields that we are able to build in the funnel at these 
resolutions and the short evolution times. Second, the zero- 
conductivity plasma in the funnel is just adjacent to the high- 
conductivity plasma of the torus; this large jump in the con¬ 
ductivity helps in preserving the magnetic-jet structure and in 
providing a natural agent for the collimation of the flow at low 
latitudes. 

It is useful to remark that in close analogy with what found 
in Ref. [36] the magnetic-jet structure produced here is not a 
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FIG. 6. Three-dimensional snapshots of the norm of the toroidal magnetic field |Btor| (left panel), of the poloidal one \Bpo\\ (middle panel), 
and of the total magnetic field \B\ (right panel), all at time t = 18.3 ms. Additionally, we plot the modulus of the magnetic field at the same 
time to illustrate where the toroidal and poloidal contributions become dominant. Note that the magnetic field at the edges of the funnel starts 
developing a toroidal component that exhibits signs of twisting, while the magnetic field in the evacuated region is predominantly poloidal. 
The figures show a cut through the torus, with the z-axis facing down, in order to be able to see the evacuated region formed along the z-axis 
and the funnel-wall structure that develops at the interface with the interstellar medium. The domain plotted corresponds to a rectangular grid 
with dimensions [0 km, 115.8 km] x [—115.8 km, 115.8 km] x [0 km, 92.16 km] 
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FIG. 7. Three-dimensional snapshots of the rest-mass density p (top left panel), specific internal energy e (top right panel), and modulus of 
the magnetic field \B\ (bottom panel) at t = 18.3 ms. Additionally, we plot the magnetic-field lines on top of these quantities to illustrate the 
topology of the magnetic field. Because of the symmetries applied in our simulation we only show a quadrant of the black hole-torus system. 
The figures show a cut through the torus, with the z-axis facing down, in order to be able to see the evacuated region formed along the z-axis 
and the funnel-wall structure that develops at the interface with the interstellar medium. The domain plotted corresponds to a rectangular grid 
with dimensions [0 km, 115.8 km] x [—115.8 km, 115.8 km] x [0 km, 92.16 km]. 


relativistic outflow. Instead, it can just be viewed as an al¬ 
most quasistationary magnetic structure confining the tenuous 
plasma in the funnel and confining it away from the dynamics 
of the ultradense plasma in the torus. In fact, despite the re¬ 
sistive losses, the plasma in the funnel does not have yet suffi¬ 
cient internal energy to be able to launch a relativistic outflow. 
It is possible that the strong magnetic flelds in the vicinity of 
the black hole could provide the conditions for electromag¬ 
netic extraction of the black hole’s rotational energy through 
the Blandford-Znajek mechanism [101] or through a general¬ 


ized Penrose process [102]. Alternatively, the energy required 
for launching a relativistic outflow could also be efficiently 
deposited along the baryon-poor funnel by reconnection pro¬ 
cesses not fully modelled here or by neutrino pair annihilation 
[103, 104]. Clearly, additional work is needed to assess the 
robustness of our modelling of the funnel region and to assess 
whether and how energy can be deposited in the magnetic-jet 
structure. 

A closer look at the magnetic-jet structure is offered in 
Fig. 5, which shows two-dimensional snapshots on the (x, 
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planes of the rest-mass density (top row) and of the magnetic 
field. The two columns refer to t = 11.60 ms (left column), 
when the HMNS is just about to collapse, and tot = 18.54 ms 
(right column), when the black hole-torus system is toward 
reaching a quasi stationary equilibrium. The different panels 
in Fig. 5 should be compared with the corresponding ones in 
Figs. 3 and 6 but are represented here on much larger spa¬ 
tial scales (see also Appendix A for a closer comparison of 
the magnetic-field structure in the IMHD simulations). Inter¬ 
estingly, the magnetic-jet structure extends well beyond the 
scale of direct infiuence of the black hole and shows a coher¬ 
ent structure on scales of > 250 km. The scale of the mag¬ 
netic structure is much larger than that of the accreting torus 
(i.e., ^ ±60 km) and of the black hole (i.e., ~ ±5 km). 

Additional information on the magnetic-field topology can 
be appreciated by considering three-dimensional views of the 
magnetic-field strength and field lines. This can be seen in 
Figs. 6, where we show a three-dimensional snapshot of the 
toroidal (left panel), poloidal (middle panel), and total mag¬ 
netic field (right panel). It is important to remark that a well- 
defined magnetic-jet structure has been recently reported also 
in Ref. [21] from IMHD simulations of the merger of a black 
hole-neutron-star binary. In addition to a jet structure not very 
different from the one reported here, the authors in Ref. [21] 
are also able to produce a sustained outflow from the accre¬ 
tion torus. At the same time, the IMHD simulations reported 
in Ref. [45], where a very high spatial resolution was used, 
do not reveal the formation of such a magnetic-jet structure. 
It is difficult to assess at the moment the origin of these dif¬ 
ferences, partly because of the limited amount of information 
provided on the simulations in Ref. [45]. A more extended 
discussion of the properties of the magnetic-field topology and 
dynamics, made along the lines suggested in this paper (see 
also Sec. IV C 1) would be useful to clarify if there are really 
differences and their origin. 

When studying the magnetic-field topology, i.e., whether 
it is mostly poloidal or toroidal, it is unavoidable to discuss 
where certain measurements are made, as the magnetic field 
can be at the same time mostly poloidal and mostly toroidal 
but in two different regions. A quick inspection of the three 
panels suggests that the magnetic field in the low-density fun¬ 
nel is predominantly poloidal (clearly shown in the middle 
panel of Fig. 6). It is quite natural to expect that the magnetic 
field will be essentially poloidal near the rotation axis, where 
matter has a specific angular momentum that is intrinsically 
small. Equally natural is to expect that a toroidal component 
will start to develop away from the axis and as one approaches 
the regions filled by the torus. This behavior can be explained 
by the fact that some of the magnetic-field lines in the funnel 
are anchored to the highly conducting material at the edges 
of the torus, which is rotating at nearly Keplerian velocities. 
Indeed, the left panel of Fig. 6 shows that the magnetic field 
acquires a toroidal component near the edges of the funnel, 
and then becomes essentially toroidal in the torus. This is also 
shown in Fig. 7, which offers three-dimensional snapshots of 
the rest-mass density p (left panel), of the specific internal en¬ 
ergy e (middle panel), and of the modulus of the magnetic 
field \B\ (right panel) at t = 18.3 ms. Also reported are the 


magnetic-field lines, of which the three-dimensional represen¬ 
tation confirms that the magnetic field is mostly poloidal in 
the magnetic-jet structure, acquiring a twist and a kink when 
it reaches the edges of the funnel, and becoming essentially 
toroidal inside the torus. 

Finally, we note that a vigorous outfiow develops at the in¬ 
terface between the magnetic-jet structure and the torus. The 
resulting shearing boundary layer could be the site for the 
development of a Kelvin-Helmholtz instability, which unfor¬ 
tunately we cannot investigate at the present resolutions and 
without a more sophisticated treatment of the conductivity in 
the transition between large and small values (the instability 
does potentially develop in a particularly difficult region). It is 
clear, however, that the dynamics along the torus walls has all 
the potential of yielding interesting observational features and 
should be investigated in the future, possibly using advanced 
numerical techniques such as those presented in Ref. [105]. 

A word of caution should be spent before concluding this 
section. While the behavior described above appears reason¬ 
able and is possibly the expected one on the basis of rather 
simple considerations, this behaviour is ultimately the result 
of our choice for the conductivity profile in Eq. (22) and of 
our choice of rest-mass density in the atmosphere. The large 
computational costs associated with these simulations prevent 
us from presenting at this point a systematic investigation of 
how sensitive the results are on the choice for a and Patm- We 
are aware that this represents a limitation of our investigation, 
which we plan to resolve with future simulations. 

C. Comparison with IMHD simulations 

As described in the previous section, the dynamics of the 
RMHD simulations is rather similar, at least qualitatively, to 
the one observed in IMHD simulations of the same binary pre¬ 
sented in Ref. [35]. Yet, there are some important differences, 
and these can be best appreciated if we perform a careful com¬ 
parison of the two evolutions. To this scope, we have per¬ 
formed additional simulations of the same binary discussed 
in the previous section when, however, the set of equations 
solved are those of general-relativistic IMHD. We note that 
this was necessary because the simulations in Ref. [35] used 
a different grid structure and resolution but also investigated 
quasicircular initial data in contrast to the reduced linear mo¬ 
menta we have considered here. 

In the following we present the results of this side-by-side 
comparison, first using two-dimensional spacetime diagrams 
and then moving to standard one-dimensional snapshots. 

1. Spacetime diagrams 

We have found that a very efficient way of carrying out a 
comparison between two MHD evolutions which are qualita¬ 
tively similar is to use two-dimensional spacetime diagrams. 
This technique, which was first introduced in Ref. [106], pro¬ 
vides a color-coded evolution of various scalar quantities as 
measured along principal axes (e.g., the x- and 2 ;-directions) 
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FIG. 8. Two-dimensional spacetime representation of the evolution of the rest-mass density p (left panel), of the specific internal energy e 
(middle panel), and of the magnetic field norm \B\ (right panel). All quantities are measured along the x-axis and are therefore representative 
of motions on the equatorial plane. Indicated with a solid green line is the evolution of the apparent horizon; note that the values in the color 
bars are saturated and do not correspond to the minimum and maximum values of the corresponding fields. 


and has the advantage of summarizing simply even rather 
complex dynamics. 

As representative examples, we show in Figs. 8 and 9 the 
differences between the RMHD and IMHD implementations 
in the evolution of the rest-mass density p (top left panel), 
of the specific internal energy e (top right panel) and of the 
magnetic field norm \B\ (bottom panel). Note that each panel 
is split into two diagrams, with the left one referring to the 
RMHD solution and the right one to the IMHD solution. Fur¬ 
thermore, while Fig. 8 refers to quantities measured along the 
x-axis and hence is representative of motions on the equatorial 
plane, Fig. 9 refers to the z-axis and is therefore representa¬ 
tive of motions in the polar region. In all panels we indicate 
with a solid green line the evolution of the apparent horizon 
(see the discussion in Ref. [89] on how to interpret such a 
line). Note that the values in the color bars are saturated and 
do not correspond to the minimum and maximum values of 
the corresponding fields. 

It is clear from both figures that the evolution of all quanti¬ 
ties is very similar during the inspiral, when indeed the IMHD 
and RMHD evolutions should be mathematically identical 


given that our fields are contained in the stars. As the binary 
reaches the merger, matter is expelled from the stars mainly 
in the equatorial plane, as is evident in Fig. 8, which also 
shows that the ejected matter moving through the low-density 
medium generates shocks that heat up the plasma and appear 
as thick lines. The subsequent bursts of matter are mainly 
associated with the fundamental mode (/-mode) of oscilla¬ 
tion of the HMNS [13, 107]. The equatorial ejections are also 
accompanied by four or five subsequent “bursts” along the 2 ;- 
axis, starting in both simulations at approximately 5 ms, and 
can be seen in Fig. 9. Clearly, the violent oscillations experi¬ 
enced by the HMNS launch matter essentially isotropically. 

We have already mentioned that the differentially rotating 
HMNS does collapse promptly to a black hole by first re¬ 
arranging its angular velocity profile and mass distribution 
through magnetic braking. During this phase, angular mo¬ 
mentum is transported outward, with the outer fluid elements 
moving further away from the star and the inner ones mov¬ 
ing toward the center as a result of the angular-momentum 
losses. Clearly, this redistribution of angular momentum will 
be different in the RMHD and IMHD evolution, with the latter 
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FIG. 9. The same as in Fig. 8 but for quantities measured along the z-axis and therefore representative of motion in the polar regions. 


being more efficient in transporting angular momentum out¬ 
ward (in IMHD the fluid can only move along magnetic-field 
lines, while it can also partially cross them in RMHD). As a 
result, the collapse can take place slightly earlier, occurring at 
11.9 ms in the IMHD evolution and at 13.8 ms in the RMHD 
simulation (cf.. Figs. 8 and 9). 

We have already commented that one of the major differ¬ 
ences between the IMHD and RMHD simulations is that our 
implementation of the latter allows for a description of prop¬ 
agating electromagnetic waves in vacuum. By contrast, the 
atmosphere treatment of the IMHD implementation is such 
that it does not evolve the magnetic field as the fiuid veloci¬ 
ties are reset to zero there. As a result, already after the first 
0.5 ms the magnetic field manages to diffuse out of the stars, 
heating up the plasma in the outer layers and forcing it to ex¬ 
pand. In this way, magnetic energy is converted into internal 
energy, and therefore the magnetic field at the center of the 
stars is (slightly) lower than in the corresponding IMHD sim¬ 
ulation. In contrast, the magnetic field along the z-axis in the 
RMHD simulation is higher than in the IMHD simulation in 
the first few milliseconds. Indeed, the differences between the 
IMHD and RMHD simulations become particularly evident 
after black-hole formation, when the funnel is evacuated and 
the magnetic-jet structure is built (see the bottom panels of 


Figs. 8 and 9). The magnetic diffusivity in the RMHD simu¬ 
lation acts so rapidly that in a bit more than one crossing time 
the whole computational domain in the RMHD run is filled 
with electromagnetic fields despite the fact that the magnetic 
field was initially constrained to the stellar interior.^ 

Note that the modulus of the magnetic field along the 2 ;- 
axis is about 2 orders of magnitude larger in the RMHD sim¬ 
ulation (cf., - 10^^ G in IMHD vs - 10^^ G in RMHD). 
This is due (mainly) to the intense currents produced by the 
rapidly rotating torus and (partly) to the magnetic-field diffu¬ 
sion of the strong magnetic field in the torus, which diffuses 
across the walls of the funnel. This behavior of the magnetic 
field is in fact similar to the one observed in a stable mag¬ 
netized star with extended magnetic fields which was stud¬ 
ied in Ref. [66], where after an initial transient, the system 
relaxed to a solution consisting of a large-scale, nearly elec¬ 
trovacuum, dipolar magnetic field configuration in the exte¬ 
rior, which was anchored to the highly conducting neutron 


^ We prefer to be repetitive rather than confusing: the magnetic field diffuses 
out of the stellar matter because of the finite resistivity at the surface of the 
neutron stars, of the HMNS, or of the torus. Once in the atmosphere, how¬ 
ever, the magnetic fields propagate as electromagnetic waves in vacuum. 
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FIG. 10. Evolution of the maximum of the rest-mass density for 
both the IMHD (blue solid line) and the RMHD simulation (green 
solid line). The top panel refers to the initial data with reduced linear 
momenta, while the bottom one refers to data in quasicircular orbits. 
Also shown in the top panel is the evolution of the RMHD set of 
equations with a large and uniform conductivity (red solid line); in 
this case the evolution should mimic the IMHD one, as indeed it 
does. 


star. The larger values of the magnetic field in the funnel are 
particularly encouraging as strong magnetic fields are neces¬ 
sary to produce a large acceleration along the 2 ;-direction. The 
maximum Lorentz factors achieved after the collapse in both 
simulations is TF « 2.0 — 2.7, with the highest values occur¬ 
ring at the end of the simulations. 

The magnetic field in the torus is also stronger in the 
RMHD simulation than the corresponding field in the IMHD 
run, although the differences in this case are only of 1 order 
of magnitude. The reason for this has to be found in the fact 
that in the RMHD simulation the torus is more massive and 
therefore able to sustain larger amounts of magnetic fields (as 
for isolated stars, also a self-gravitating torus in MHD equi¬ 
librium will be able to sustain stronger magnetic fields for in¬ 
creasing masses; cf., Table I). 


2. Angular-momentum transfer and HMNS lifetime 

The top panel of Fig. 10 reports the evolution of the maxi¬ 
mum of the rest-mass density Pmax for both the IMHD (blue 
solid line) and the RMHD simulation (green solid line). After 
the merger takes place at f 3.5 ms, the maximum rest-mass 
density oscillates at the /-mode frequency and experiences a 
sudden drop when an apparent horizon is found since matter 
inside the apparent horizon is excluded from the calculation of 
Pmax- It is quite obvious that the behavior of Pmax is different 


in the two simulations, with the HMNS in the RMHD surviv¬ 
ing for a longer time. Note, however, that the first increase 
in the central rest-mass density in the top panel of Fig. 10 is 
slightly larger in the RMHD case, most likely because of the 
resistive increase of the internal energy at the merger. 

Also shown in the top panel of Fig. 10 is the evolution of 
a simulation in which the set of RMHD equations is used to¬ 
gether with a large and uniform conductivity a = ao = 10^ 
(red solid line). In this case, despite the very different set of 
equations solved and the different approach for enforcing the 
divergence-free condition of the magnetic field, the RMHD 
evolution should mimic the IMHD one. This is indeed the 
case, and it provides considerable confidence on the robust¬ 
ness of our RMHD approach. It suggests that the delayed 
collapse is an effect associated with the choice of physical 
resistivity and not due to numerical artifacts. 

Despite the complex dynamics, the differences between the 
IMHD and RMHD runs are not difficult to explain. As men¬ 
tioned in the previous section, in fact, an important difference 
between the IMHD and RMHD simulations is that in the latter 
the magnetic field cannot be perfectly locked with the plasma. 
As a result, the IMHD evolution is more efficient in redis¬ 
tributing the angular momentum in the system and, in particu¬ 
lar, in transporting it outward. This magnetic-braking process 
deprives the HMNS core of the angular-momentum support, 
and this leads to an earlier collapse. Clearly, since the conduc¬ 
tivity in the HMNS interior is very high also in the RMHD 
simulation (although not infinite), magnetic fiux freezing is 
very efficient here as well and the differences in the dynam¬ 
ics of the IMHD and RMHD simulations can only be small. 
This explains why overall the time of collapse varies by only 
^ 1.9 ms. In addition, the important differences in this dy¬ 
namics are also expected to take place in the outer layers of 
the HMNS, where the conductivity decreases as a response 
to the conductivity profile (22). We should note that the dif¬ 
ference in the lifetime of the HMNS is not related to the use 
of initial data with modified momenta, but is present also for 
a binary of which the initial data is on a quasicircular orbit. 
This is shown in the bottom panel of Fig. 10, where the corre¬ 
sponding quantities are shown, and where it is clear that also 
in this case the HMNS collapses earlier to a black hole in the 
IMHD simulation. We also recall that for the binary in qua¬ 
sicircular orbit no additional refinement level is added after 
the merger. Hence, when comparing in Fig. 10 the RMHD 
evolution of the binaries with reduced momenta and on quasi¬ 
circular orbits, one is also effectively comparing the evolution 
of the HMNS at two different resolutions, i.e., 148 and 296 m, 
respectively. The fact that they both yield a longer lifetime of 
the HMNS provides an indirect validation of the numerical 
consistency of the RMHD solution. 

Because the differences in the magnetic braking between 
the RMHD and IMHD implementations are intrinsically 
small, it is not easy to show that it is exactly these differences 
that are responsible for the earlier collapse of the HMNS in 
the IMHD simulations. However, such evidence is offered 
in Fig. 11, which reports the spacetime diagrams along the 
x-direction of the specific angular momentum i := —U(^lut 
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FIG. 11. Spacetime diagrams of the specific angular momentum i := —Uc^/ut for both the RMHD (left panel) and IMHD (middle panel) 
simulations. In addition, we show the relative difference in the specific angular momenta between the RMHD and IMHD implementations 
(right panel). 


[68].^^ The left panel refers to £r, the specific angular mo¬ 
mentum of the RMHD run, while the middle panel refers to 
the corresponding quantity for the IMHD simulation, ij, and 
the right panel refers to the relative difference: 1 — IrIIi. 

A rapid inspection of the left and middle panel reveals that 
in both the RMHD and IMHD simulations the specific angu¬ 
lar momentum increases outward (as it should for a rotating 
fluid satisfying the Rayleigh stability criterion), but also that 
the profiles are not constant in time and show instead periodic 
variations that are in phase in the two simulations. These vari¬ 
ations reflect the large oscillations of the HMNS, and, indeed, 
the oscillations in Ir, ij take place at the same frequency as 
those in the rest-mass density and shown in the other space- 
time diagrams in Fig. 8 (note the different scale in the x-axis). 
However, there are important small differences in the dynam¬ 
ics of Ir and ij, which are apparent in the right panel of 
Fig. 11, where regions in red indicate that the specific angular 
momentum in these regions is higher (of ^ 20 — 30%) in the 
IMHD simulation than in the RMHD simulation, while blue 
regions are exactly the opposite. It is apparent that the differ¬ 
ences are larger in the central parts of the HMNS, while the 
specific angular momenta are very similar in the outer layers, 
i.e., for X > 5 km. We recall that an excess of specific angular 
momentum at a given position on the x-axis reflects fluid ele¬ 
ments that are rotating at a larger frequency (i = VLx‘^ in the 
Newtonian limit) and this is indeed what one would expect if 
the magnetic fields and the fluid are tightly coupled. Hence, 


A very similar behaviour is observed if the specific angular momentum is 
shown in terms of ^ hU(j) (see [68] for a discussion in the differences in 
the two definitions). 


the red regions in the right panel Fig. 11 can be taken to signal 
a more efficient transfer of angular momentum from the inner 
regions of the HMNS. A direct consequence of this transfer 
of angular momentum is the appearance of blue regions ad¬ 
jacent to the red ones and signalling therefore fluid elements 
that have slowed down. 

Although the differences in £r and £i are small, the trans¬ 
fer continues steadily and with an increased rate up to f « 
11 ms, when a much larger transfer of angular momentum 
takes place. This signals the onset of the instability to grav¬ 
itational collapse in the IMHD simulation, which effectively 
takes place soon after, i.e., at t 12 ms. 

Three remarks should be made before concluding this sec¬ 
tion. First, the fact that an RMHD simulation with a uniform 
conductivity yields the same collapse time as an IMHD sim¬ 
ulation gives us great confidence about the correctness of the 
RMHD evolution with nonuniform conductivity. Second, al¬ 
though the difference in the survival time between the RMHD 
and IMHD evolution is here rather small, it can be much larger 
if smaller values of the resistivity are chosen for the stellar in¬ 
terior and for less massive HMNSs; unfortunately present as¬ 
tronomical observations do not set any stringent constraint on 
the values of the conductivity at these temperatures, rest-mass 
densities, and magnetic fields. More importantly, however, a 
longer survival time is a useful new result in the modelling of 
binary neutron stars, as it points out that the HMNS can sur¬ 
vive on comparatively longer time scales than those computed 
so far in pure hydrodynamics or in IMHD simulations. This 
is not a minor detail as the most recent modelling of SGRBs 
with an extended x-ray emission invokes the existence of a 
magnetized HMNS that is able to survive on time scales of 
the order of 10^ — 10^ sec before collapsing to a black hole 

















19 


[57-62]. Finally, computational constraints have prevented us 
from extending much past black-hole formation the evolution 
of the RMHD/IMHD simulations with initial data on quasicir¬ 
cular orbits. Nevertheless, the fact that already the dynamics 
of the HMNS is unaffected by the initial reduction in linear 
momenta and that the HMNS collapses to a black hole earlier 
in both cases, provides us with confidence that the magnetic 
field dynamics discussed in Sec. IV B will be very similar also 
for binaries on quasicircular orbits. This will also be the focus 
of our future work. 


5. magnetic field growth 

In Secs. IV B and IV C 1, we have already discussed the 
properties of the evolution of the magnetic fields, but have not 
quantified in detail how the magnitude of the magnetic field 
changes with time and how this evolution varies in the RMHD 
and IMHD simulations. This is done now in Fig. 12, where 
we report the evolution of the maximum of the modulus of the 
magnetic field (black lines) but also of its toroidal (red lines) 
and poloidal (blue lines) components, either in the RMHD 
simulation (left panel) or in the IMHD simulations (middle 
panel). 

Note that the evolution of the magnetic field in the two 
implementations is very similar during the inspiral, but also 
that this changes considerably after the merger. While in 
both cases the toroidal magnetic field grows exponentially, the 
growth is of about 1 order of magnitude in the IMHD simula¬ 
tion but of a factor 2 smaller for the RMHD simulation. Fur¬ 
thermore, the magnetic field reaches values of 3 x 10^^ G just 
before the collapse in the IMHD simulation and 8.3xlO^^G 
only just after the merger in the RMHD simulation. This dif¬ 
ferent behavior is not difficult to explain and is simply due to 
the fact that the shearing of magnetic-field lines is less effi¬ 
cient in RMHD because of the finite conductivity of the mat¬ 
ter. Because the growth at the merger is mostly due to the 
shear layer between the two impacting stars, it is quite nat¬ 
ural that a resistive calculation will lead to a smaller mag¬ 
netic field, quite independently of how well the instability is 
resolved. 

As for the rest-mass density (cf.. Fig. 10), the collapse of 
the HMNS to a black hole leads to a rapid decrease of the max¬ 
imum value of the magnetic field, as shown in Fig. 12 where 
the vertical lines signal the first appearance of the apparent 
horizon. Also in this case, the strongest magnetic fields are 
hidden inside the horizon, and the maximum values reported 
are those relative to the magnetic field in the torus. As re¬ 
marked already when commenting on the spacetime diagrams 
in Sec. IV C 1, the larger values of the magnetic field in the 
RMHD simulation are the result of a more massive torus pro¬ 
duced in this case (cf.. Table I). 

In the right panel of Fig. 12 we complement the evolution 
of the magnetic fields with the evolution of the I/ 2 -norm of 
the divergence of the magnetic field for the RMHD simula¬ 
tion (green dashed line) and for the IMHD one (blue solid 
line) after it is multiplied by 10^^. We recall that the IMHD 
simulation makes use of a constrained transport scheme [108], 


and hence it is able to maintain the violations of this con¬ 
straint down to machine precision. The RMHD simulation, 
on the other hand, makes use of a divergence-cleaning scheme 
[67], which is widely known to be less efficient in suppress¬ 
ing the violations. Yet, the purpose of making this com¬ 
parison is mostly that of highlighting that the divergence¬ 
cleaning approach used here may not be as efficient as the 
constrained transport, but it yields nevertheless very small 
violations. Indeed, the evolution of the I/ 2 -norm of the ra¬ 
tio between the divergence and the magnetic field strength, 
i.e., \\ViB^/^/BiB '^\\2 and not shown in Fig. 12, is of the 
order of ^ 10“^, thus similar to the values reported by 
similar works in the literature [87]. Note that the late-time 
moderate growth of the divergence in the RMHD simula¬ 
tion is due to the amplification of the magnetic fields in the 
torus, and it would be interesting to investigate if a dynami¬ 
cally adapted dissipation parameter for the divergence clean¬ 
ing method could help reduce such a growth. 


4. Magnetically driven wind and bursting activity 

In Secs. IVA and IVCl, we have already anticipated 
that a wind is produced after the merger, either as a result 
of shock heating at the merger or because of magnetic wind¬ 
ing and consequent pressure imbalance in the outer layers of 
the HMNS or because of neutrino losses [109, 110]. In ad¬ 
dition to an almost quasistationary and quasi-isotropic wind, 
both the RMHD and the IMHD simulation show the existence 
of mildly anisotropic and quasiperiodic launching of low rest- 
mass density, high internal energy blobs of matter that we will 
refer to as bursts. Overall, seven bursts are launched during 
the total time of the simulation, with five bursts relative to the 
HMNS stage and two being produced after black-hole forma¬ 
tion. 

In particular, the first two bursts eject material that is mov¬ 
ing through the low-density atmosphere with an average speed 
of ^ 0.4 — 0.6 but that decelerate as they move away from the 
HMNS, reaching a final outward velocity of ^ 0.16 — 0.18. 
This behavior is in part due to the natural conversion of ki¬ 
netic energy to binding energy but also to the interaction of 
the bursting material with the slow isotropic wind. This inter¬ 
action, which is obviously accompanied by shocks, provides a 
damping mechanism on the propagation of the ejected mate¬ 
rial. However, the slow down of the baryon rich material could 
be revived later on if the slow wind is impacted by a faster, 
baryon poor wind, as suggested in the “two-winds” model of 
Ref. [61]. 

We have investigated the properties of the bursts by study¬ 
ing the evolution of the specific internal energy as measured 
by an observer on the z-axis. This is reported in Fig. 13 for an 
observer at (x, y, z) = (0, 0, 86.82) km, for both the RMHD 
simulation (green dashed line) and the IMHD one (blue solid 
line). Clearly, after the merger the specific internal energy 
exhibits a quasiperiodic behavior in both simulations, with a 
very rapid increase. The increase in the specific internal en¬ 
ergy is of slightly less than 2 orders of magnitude and is fol¬ 
lowed by a slower decay. 
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FIG. 12. Evolution of the maximum of the magnetic-field strength |B|max for both the RMHD (left panel) and the IMHD simulation (middle 
panel). The black solid lines correspond to the time series of the maximum magnetic-field modulus; the blue dashed lines correspond to 
the time series of the maximum poloidal magnetic-field norm, |B|max; and the red dotted-dashed lines correspond to the maximum toroidal 
magnetic field, |max. Additionally, we show the L 2 -norm of the divergence of the magnetic field, for the RMHD simulation (green dashed 
line) and IMHD one (blue solid line). The collapse times are depicted with black dotted or black dashed lines. 



FIG. 13. The time series of the specific internal energy measured at 
(x, y, z) = (0, 0, 86.82) km are shown for the RMHD simulation 
with a green dashed line and for the IMHD simulation with a solid 
blue line. The time of collapse for the RMHD case is depicted with 
a black dashed line, while it is depicted with a black dotted line for 
the IMHD one. 


As already mentioned in Sec. IVCl, this behavior can 
be associated with the oscillations of the HMNS and is ob¬ 
served also in the evolution of the rest-mass density (cf., the 
five peaks in Fig 10 and in Fig. 9, top left panel). The char¬ 
acteristic frequency of these peaks is related to the /-mode 
frequencies of the bar-deformed HMNS [13, 107] and thus 
these bursts occur every ^ 1.7 — 2.0 ms. Similar peaks (al¬ 
though less marked) can be observed also in the magnetic field 
along the 2 ;-axis (cf., Fig. 9, bottom panel). The difference is 
that they occur slightly later than the specific internal energy 
ones, possibly indicating a conversion of kinetic energy into 


magnetic energy and vice versa. The rest-mass density in the 
blobs of low-density material ejected in the bursts depends on 
height, but, at a distance ofz = 86.82 km along the z-axis, it 
is 6 X10^ gr cm“^—6 X10^ gr cm“^, while the magnetic field, 
at the same location, has a strength of^ 10^ — 5 x lO^^G. 

The bursting activity continues also after the collapse, but 
with somewhat different properties. First, the frequency is 
now set by the radial epicyclic frequencies of the oscillating 
torus as deduced, for instance, when analyzing the time series 
of the specific internal energy at a ~ 60 km on the x-axis (see 
Refs. [106, 11 1] for an introduction to these frequencies in ro¬ 
tating tori). Second, the rest-mass density of the blobs ejected 
is smaller and of the order of 5 x 10^ gr cm“^, while the mag¬ 
netic field oscillates between 10^ —10^ G and is stronger prob¬ 
ably as a result of the magnetic field increase in the funnel. 

With only two bursts observed after black-hole formation, 
the time span is too short to reach a firmer conclusion on 
the origin of the bursts after the collapse. However, all 
present evidence seems to suggest that the postmerger bursts 
are triggered by an increased mass-accretion rate as the torus 
approaches the black hole during the inward phase of its 
epicyclic oscillation. Of course also resistive reconnection 
processes could lead to the conversion of magnetic energy into 
internal energy and thus may be invoked to explain this phe¬ 
nomenology. We find this not a likely explanation mostly be¬ 
cause of the ordered magnetic field structure that builds up in 
the funnel and which seems rather stationary. Clearly, a more 
detailed study of long-term evolutions with different prescrip¬ 
tions for the conductivity profiles is necessary before reaching 
more robust conclusions. 

We finally note that the outflows produced either by shock¬ 
heating, magnetically driven winds or by the periodic bursts, 
eject a substantial amount of matter. More speciflcally, the to¬ 
tal rest-mass flux across a spherical surface located at r = 
295.4km is measured to be 0.5 — 2.0 Mq sec“^, which 
amounts to a total of ^ 0.01 Mq ejected from the beginning 
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M (Mo) J/M^ Mtor. (Mo) 

rtor. (km) 

RMHD 

2.88 0.873 0.095 

105.9 

IMHD 

2.91 0.884 0.075 

88.9 


TABLE I. Properties of the black hole-torus system at t = 4.74 ms 
after the appearance of the apparent horizon for both the resistive 
and IMHD simulations. Shown are the mass and dimensionless spin 
of the black hole, as well as the rest mass and size of the torus as 
estimated with a cutoff on the rest-mass density at p = 10^°g cm“^. 


of the simulation and over the survival time of the HMNS 
(i.e., 13.8 ms). We also note that the mass-ejection rates re¬ 
ported here are larger than those obtained in Ref. [43], where 
mass fluxes of^ 10“^ — 10“^ Mq sec“^ were reported. 
However, this is not particularly surprising and for a num¬ 
ber of reasons. First, the HMNS considered here is the self- 
consistent result of a binary merger, while the one studied in 
Ref. [43] was built using an axisymmetric differentially rotat¬ 
ing equilibrium with a standard (but somewhat arbitrary) law 
of differential rotation. Second, as mentioned in Sec. HID, 
the linear momenta in our initial data are artificially reduced 
to accelerate the inspiral. This also leads to a more violent 
merger and to larger mass losses at least till black-hole for¬ 
mation. After the HMNS collapses, in fact, the mass fiux sat¬ 
urates at 0.2 Mq sec“^. Finally, the initial magnetic field in 
Ref. [43] is about 2 orders of magnitude larger, and this facili¬ 
tates substantially the loss of MHD equilibrium at the surface 
of the HMNS and thus the mass loss. 


5. Black hole-torus properties 

Because the magnetic field is not strong enough to alter the 
torus dynamics significantly [35], it is natural to expect the 
two solutions (IMHD and RMHD) to be very similar in terms 
of dynamical properties once a nearly quasistationary state is 
established after the collapse to a black hole. This expectation 
is confirmed in the upper panel of Fig. 14, which shows the 
rest-mass density distribution along the x-direction at about 
4.74 ms after the formation of the apparent horizon. Clearly, 
the two distributions are very similar but not identical. The 
RMHD simulation, in particular, yields larger rest-mass den¬ 
sities in the outer portions of the torus, which, in turn, are 
responsible for larger rest masses (cf., Table I) and stronger 
magnetic fields (see the discussion in Sec. IV C 1). Further¬ 
more, the RMHD simulation also yields a larger torus as mea¬ 
sured from the average position on the x-axis where the rest- 
mass density falls below p = lO^^g cm“^. 

Additionally, the bottom panel of Fig. 14 illustrates the an¬ 
gular velocity profile along the x direction ft := u^/u^ for 
the RMHD implementations (green dashed line) and for the 
IMHD one (blue solid line), at the same time as the upper 
panel. Also shown as a black dashed line is a reference Keple- 
rian profile, i.e., scaling like and which is well matched 
by both distributions. As remarked in Ref. [106], the fact that 
the outer parts of the torus have a quasi-Keplerian behavior 



FIG. 14. Upper panel, rest-mass density profile along the x-axis 
for the RMHD (green dashed line) and IMHD (blue solid line) sim¬ 
ulations at time t = 4.74 ms after the apparent-horizon formation. 
Bottom panel the same as above but for the angular velocity; also 
shown as a reference is a Keplerian profile (black dashed line). 


has two important implications. First, it suggests that the tori 
will be stable and not subject to dynamical instabilities that 
would lead to their rapid destruction (see, e.g., Refs. [112- 
114]). Second, a quasi-Keplerian profile also provides optimal 
conditions for the development of an MRI in the torus, thus 
opening the possibility of further amplification of the mag¬ 
netic fields that are present after the collapse of the HMNS. 


6. Electromagnetic luminosities 

Despite the exploratory nature of the simulations carried 
out here, we have computed for both the RMHD and IMHD 
implementations the total electromagnetic luminosity 
emitted. This has been estimated as a surface integral of the 
Poynting fiux over spherical surfaces placed at representative 
coordinate radii r^, where has been varied to guarantee 
that the measurement is an asymptotic one and is not affected 
by the local plasma dynamics. We recall, in fact, that, because 
in the IMHD approximation the magnetic fields are locked 
with the plasma, the electromagnetic luminosity estimates can 
be heavily influenced by the presence of matter in the outer 
regions and thus not correspond to a genuine amount of elec¬ 
tromagnetic energy fiux leaving the system. Unfortunately, 
there is no simple way within the IMHD approximation of de¬ 
termining whether the integral of the Poynting fiux computed 
on the numerical grid is genuinely asymptotic. However, it is 
certainly reassuring if the values of are independent from 
the extraction radius. 

The evolution of the electromagnetic luminosity is illus¬ 
trated in Fig. 15, where we report it as computed at three dif¬ 
ferent extraction radii, = {147.71, 221.57, 295.43} km. 
Note that the luminosity during the inspiral phase of the 
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FIG. 15. Time series of the Poynting flux computed at different extraction radii = {147.71, 221.57, 295.43} km are shown here for the 
RMHD implementation (left panel) and the IMHD one (right panel). The dotted line represents the time of the merger, and the black solid line 
represents the time of horizon formation. 


RMHD simulation is much larger than the corresponding one 
in the IMHD simulation because of the diffusion of the mag¬ 
netic field across the stars’ surfaces, that fills the entire do¬ 
main with vacuum electromagnetic fields. As remarked al¬ 
ready in Sec. IIIC 1 and IV C 1, these magnetic fields are in 
areas which are treated as atmosphere from a hydrodynami- 
cal point of view, but where the Maxwell equations are solved 
with zero conductivity, so that the electromagnetic fields can 
propagate freely. 

We note that a nominal value of a = 0 does not imply that 
the electromagnetic luminosity will be zero, since the motion 
of the compact stars will introduce perturbations in the ex¬ 
ternal electric and magnetic fields, and thus a net Poynting 
flux (see Ref. [115] for the electromagnetic emission of inspi- 
ralling binary black holes in electro vacuum). Indeed, we have 
found that the electro vacuum luminosity before the merger is 
^ sec“^, which is smaller than the one reported 

in Ref.[56] (i.e., ^ 10^^ erg sec“^), where the stellar 

exteriors were modelled in the force-free approximation. Al¬ 
though it has already been found that the electrovacuum lumi¬ 
nosity is slightly smaller than the force-free one for the same 
system (see Refs. [116-118]), the differences found here are 
larger than expected and this may be due to the rather differ¬ 
ent way in which the exterior regions of stars are treated. By 
contrast, the electromagnetic luminosity before the merger in 
the IMHD simulation, where the magnetic fields are always 
contained inside the stars, is essentially zero. 

After the merger, the electromagnetic luminosity grows 
rapidly of about 2 orders of magnitude, essentially as a result 
of the growth of the magnetic field already discussed in the left 
panel of Fig. 12 (we recall that the electromagnetic luminos¬ 
ity should scale quadratically with the magnetic field). During 
the postmerger phases, the luminosity ranges from rsj 10^® to 


lO^i gj.g sec“^, to reach values up to 10^^ erg sec“^ after the 
collapse of the HMNS to a black hole. 

In the left panel of Fig. 15, the luminosity computed on 
a surface of radius = 147.71 km (black dot-dashed line) 
does not overlap with those computed on larger radii (dark- 
blue solid line and light-blue dashed solid lines), signalling 
that this radius is too close to the central object and contami¬ 
nated by the presence of matter. Fortunately, however, the lu¬ 
minosities at = 221.57km and = 295.43 km are very 
close to each other, confirming the robustness of these mea¬ 
surements. By contrast, the three luminosities in the IMHD 
simulation reported in the right panel of Fig. 15 provide three 
different values for the luminosity, indicating that at least two 
of them (i.e., those at the smaller extraction radii) are proba¬ 
bly contaminated by the presence of bound matter and hence 
not reasonable. 


V. CONCLUSIONS 

We have presented general-relativistic simulations of the in¬ 
spiral and merger of binary neutron stars when evolved solv¬ 
ing the coupled set of the Einstein equations and those of 
RMHD. Our main interest here has been to assess the im¬ 
pact that resistive effects have on the dynamics of these bina¬ 
ries, and which are usually investigated in the more idealized 
framework of IMHD. 

Because the differences with an IMHD description could 
be rather small in certain stages of the process (e.g., during 
the inspiral), we have carried out a close comparison between 
two simulations evolving the same binary, either in the con¬ 
text of RMHD or in that of IMHD. More specifically, we have 
studied the dynamics of an equal-mass binary system of neu- 
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tron stars with a total ADM mass =3.25 Mq and an 

initial orbital separation of 45 km. The stars are initially ir- 
rotational and with zero magnetic field. A dipolar magnetic 
field is therefore added before the evolution is started, which 
is entirely contained inside the stars, at least initially. Fur¬ 
thermore, to reduce computational costs and “accelerate” the 
inspiral, we have slightly reduced the linear momenta of the 
initial data as done in Ref. [95], so that the merger occurs after 
approximately one orbit. 

A crucial goal of our RMHD approach has been that of at¬ 
taining a smooth resistive description from the highly con¬ 
ducting, high-density stellar interior, out to regions of very 
low-density plasma, where the electromagnetic fields decou¬ 
ple from the fluid. Falling between these two regimes is the 
large amount of high-density, small-velocity material that is 
ejected during and after the merger by the HMNS and, once 
the latter collapses to a black hole, by the accreting torus. This 
material, occupies a large portion of the computational do¬ 
main and is produced either by the spiral arms launched at the 
merger, or by the magnetic winding and consequent pressure 
imbalance in the outer layers of the HMNS. Neutrino losses 
can also be a source of a wind, but we do not model this here. 

While there are several ways of potentially reaching a 
smooth transition between the IMHD limit in the stellar in¬ 
terior and an electrovacuum behavior, we have here adopted 
the same approach we have extensively investigated with iso¬ 
lated neutron stars in Ref. [66]. In essence, this matching 
is achieved through a carefully chosen conductivity profile, 
where the conductivity is directly related to the conserved 
rest-mass density and is set to zero once the latter reaches a 
value close to the atmospheric floor. This prescription has at 
least two free parameters. First, they ensure that the transition 
region covers only a thin layer close to the surface of the star. 
Second, they guarantee that this layer remains “thin” even in 
the first steps of the evolution, when the outer layers of the star 
expand due to a nonzero pressure in the atmosphere. While we 
have set these two parameters to sensible values, their infiu- 
ence on the results still needs to be fully explored. 

Overall, we have found that there are many similarities be¬ 
tween the RMHD and IMHD evolutions, but also one impor¬ 
tant difference, namely, that the survival time of the hypermas- 
sive neutron star, which increases in a RMHD simulation. The 
increased lifetime of the HMNS appears to be due to a less ef¬ 
ficient magnetic-braking mechanism in the resistive regime, 
in which matter can move across magnetic-field lines, so that 
the outward transport of angular momentum is reduced. This 
interpretation is supported by the analysis of the evolution of 
the specific angular momentum, and it shows that the transport 
is more efficient in the IMHD simulation. An extended life¬ 
time of the HMNS could have intriguing astrophysical conse¬ 
quences, since a longer-lived magnetized hypermassive neu¬ 
tron star brings support to the recent modelling of SGRBs 
in terms of long-lived magnetarlike objects produced by the 
merger [57, 58, 60, 61]. 

Another important result of these simulations is the con¬ 
firmation that a magnetic-jet structure is formed in the low- 
density funnel produced by the black hole-torus system. We 
note that these simulations have been carried out at higher res¬ 


olutions and with a different grid structure than those in Ref. 
[36]. In the RMHD simulations the magnetic-jet structure is 
far more regular, essentially axisymmetric, and extending out 
to the largest scale in our system. This is most likely the result 
of the effective decoupling established between the dynamics 
of the plasma and that of the electromagnetic fields. In the 
IMHD simulations, a magnetic-jet structure is still present, 
but on the scale of the torus. This difference is due to the 
fact that a decoupling of the electromagnetic fields from the 
plasma is not possible in the IMHD approximation, and the 
magnetic field follows tightly the turbulent dynamics of the 
matter. In this case, the magnetic-field lines are almost paral¬ 
lel to the z-axis (in analogy with what was shown in Ref. [36]) 
and the topology becomes more turbulent on large scales. In 
both regimes, the magnetic field is predominantly toroidal in 
the highly conducting torus and predominantly poloidal in the 
nearly evacuated funnel, although in the IMHD simulation, 
this happens near the rotation axis. The matter in the funnel 
does not have an internal energy sufficiently large to launch 
a relativistic outflow. However, it is reasonable to expect that 
reconnection processes or neutrino annihilation occurring in 
the funnel, none of which we model here, could potentially 
increase the internal energy in the funnel. 

The final comment of this work is in fact a caveat. While 
the dynamics of the magnetic-field results presented here ap¬ 
pears reasonable, matching the expectations for this type of 
system as well as previous simulations, we should remark that 
our results are ultimately dependent on the choice made for 
Ohm’s law and for the conductivity profile. Again, while our 
choice is a very conservative and a plausible one, it represents 
a choice nevertheless. The large computational costs associ¬ 
ated with these simulations have prevented us from presenting 
a systematic investigation of how sensitive the results are on 
the choices for Ohm’s law, for the conductivity profile, or on 
the treatment of the atmosphere. All of these issues deserve 
further investigation and will be the focus of our future work. 
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Appendix A: Magnetic-jet structure in the IMHD simulations 


Although the IMHD simulations are not the focus of this 
paper, for completeness we provide in this Appendix a rapid 
overview of the properties of the magnetic-jet structure as ob¬ 
tained within this approximation. The essence of the results 
is shown in Fig. 16, which represents the equivalent of Fig. 
5 but within IMHD. The different panels show large-scale, 
two-dimensional snapshots on the (x, z) planes of the rest- 
mass density (top row), and of the magnetic field (bottom 
row). The left column of Fig. 16 refers to t = 10.25 ms (left 
column), when the HMNS has not yet collapsed to a black 
hole, while the right column refers to t = 18.89 ms, when a 
black hole has already been formed. Because in the IMHD 
approximation the magnetic fields are tightly locked with the 
matter, it does not come as a surprise that no ordered mag¬ 
netic field structure seems to develop before the HMNS col¬ 
lapses and forms a black hole. This is because the dynamics 


of the plasma is quite turbulent at the merger and during the 
HMNS stage. However, after a black hole is formed, a well- 
ordered magnetic-field structure appears as the system reaches 
a quasistationary state. We note again that the formation of a 
magnetic-jet structure occurs around the black-hole rotation 
axis. 

Differently from the corresponding RMHD simulation, the 
magnetic-jet structure here is not very regular on large scales, 
and it is necessary to go down to the length scale of the torus, 
as shown in Fig. 17, for the magnetic-jet structure to become 
evident. Note that the magnetic-field lines are almost parallel 
to the z-axis, in analogy with what was shown in Ref. [36]. 
Finally, although we are here using only the projection of the 
magnetic-field lines on the (x, z) plane, the magnified view 
in Fig. 17 reveals that the magnetic field in the low-density 
funnel is still predominantly poloidal, although not as ordered 
as in the RMHD simulation (cf.. Fig. 3, left and right panels). 
At the same time, and not shown here for compactness, the 
magnetic field is essentially toroidal in the torus. 
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FIG. 16. The same as Fig. 5, but for the IMHD simulation. The two columns refer to t = 10.25 ms (left column), when the HMNS has not 
yet collapsed, and to t = 18.92 ms (right column), when a black hole has already been formed. Note again the formation of a magnetic-jet 
structure around the black-hole rotation axis, which however is far less regular on large scales than the one produced in the RMHD simulation. 
See also Fig. 17 for a view on smaller scales. 
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FIG. 17. The same as Fig. 16, but on a scale of [—60,60] km on the x-axis and of [0,80] km on the z-axis. Note that the magnetic-jet structure 
becomes more evident on these scales. See also Fig. 3 for the corresponding quantities in the RMHD simulation. 
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